/* Copyright (C) (2010) (Can YANG) <eeyang@ust.hk>

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public License
as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.

To collect contigency tables,
We first pre-calculate the number of 1's of 0~2^16 using the hamming weight to count the number of 1's in the string. (see http://en.wikipedia.org/wiki/Hamming_weight)
see the function bitCount. We store them in the vector "wordbits". Then use table looking method. see function popcount.

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

typedef long long   int64;
typedef unsigned long long uint64;
#define FMT_INT64   "%lld"
#define FMT_UINT64   "%llu"
#define FMT_HEX64   "%llx"

struct genotype{
	uint64 *genocase;
	uint64 *genoctrl;
};

struct MarginalDistr{
	int MarginalDistrSNP[3];
	int MarginalDistrSNP_Y[3][2];
};

static unsigned char wordbits[65536];// { bitcounts of ints between 0 and 65535 };
static int popcount( uint64 i )
{
    return( wordbits[i&0xFFFF] + wordbits[(i>>16)&0xFFFF] + wordbits[(i>>32)&0xFFFF] + wordbits[i>>48]);
}

int random(int max)
{
    int i;
    int rand_max = RAND_MAX - RAND_MAX % max - 1;
    while((i = rand()) > rand_max);
    return i % max;
}

double Abs(double a)
{
	return((a<0) ? -a : a);
}

int GetDataSize(char *filename, int **DataSize)
{
	FILE * fp, *fp_i;
	int c, ndataset;
	time_t st,ed;
	int n, p, i, flag,ii;
	char filename_i[100];


	fp = fopen(filename,"r");
	if(fp == NULL)
	{
		fprintf(stderr,"can't open input file %s\n",filename);
		exit(1);
	}

	ndataset = 0;
	while(!feof(fp)) {
		ndataset ++;
		fscanf(fp, "%s\n", &filename_i);
	}

	*DataSize = (int *)calloc( ndataset*2, sizeof(int));

	ii = 0;
	rewind(fp);
	while(!feof(fp)) {
		ii ++;
	fscanf(fp, "%s\n", &filename_i);

	fp_i = fopen(filename_i, "r");
	if(fp_i == NULL)
	{
		fprintf(stderr,"can't open input file %s\n",filename_i);
		exit(1);
	}
	printf("start getting data size of file %d: %s\n", ii, filename_i);
	time(&st);
    //initialization
	if (ii == 1)
	{
	n = 0;//samples number

	// find the number of samples: n
	while(1)
	{
		int c = fgetc(fp_i);//read a character from the data file
		switch(c)
		{
		case '\n'://the end of line
			n++;
			break;
			// fall through,
			// count the '-1' element
		case EOF://file end
			goto out;
		default:
			;
		}
	}

	}
out:
	rewind(fp_i);//Repositions the file pointer to the beginning of a file

	// find number of variables: p
	p= 0;
	i= 0;
	flag = 1;
	while(1)
	{
		c = getc(fp_i);
		if(c=='\n') goto out2;//end of line
		if(isspace(c))
		{
			flag = 1;
		}
		/*do {
			c = getc(fp);
			if(c=='\n') goto out2;//end of line
		} while(isspace(c));//space
		*/
		if (!isspace(c) && (flag==1))
		{
			p++;//indicate the dimension of the vector
			flag = 0;
		}

	}
out2:
	fclose(fp_i);

	time(&ed);

//	DataSize[0] = n;
	(*DataSize)[ndataset * 0 + ii - 1] = n;
	(*DataSize)[ndataset * 1 + ii - 1] += p-1;

	}

	fclose(fp);
	//printf("Data contains %d rows and %d column. \n", n, p);

	printf("cputime for getting data size: %d seconds.\n", (int) ed - st);
	return ndataset;
}

//see http://en.wikipedia.org/wiki/Hamming_weight

int bitCount(uint64 i)
{
	i = i - ((i >> 1) & 0x5555555555555555);
	i = (i & 0x3333333333333333) + ((i >> 2) & 0x3333333333333333);
	i = (i + (i >> 4)) & 0x0f0f0f0f0f0f0f0f;
	i = i + (i >> 8);
	i = i + (i >> 16);
	i = i + (i >> 32);
	return (int)i & 0x7f;
}


// get the marginal distribution and store it in pMarginalDistr
// pMarginalDistr[i1].MarginalDistrSNP_Y[i2][0] stores the number of cases in the cell X_i1 = i_2
// pMarginalDistr[i1].MarginalDistrSNP_Y[i2][1] stores the number of controls in the cell X_i1 = i_2
// pMarginalDistr[i1].MarginalDistrSNP[i2] stores the number in the cell X_i1 = i_2
void CalculateMarginalDistr(struct genotype *pgeno, int nsnp, int n, int nlongintcase, int nlongintctrl, struct MarginalDistr *pMarginalDistr)
{
    int i1, i2, i3;
    int count;
	//double tmp, ptmp;
    //int GenoMarginalDistr[3][2];

    for (i1 = 0; i1< nsnp; i1++)
    {
        for (i2 = 0; i2<3; i2++)
        {
            count = 0;
            for (i3 = 0; i3< nlongintcase; i3++)
			{
				count += bitCount(pgeno[i1*3 + i2].genocase[i3]);
			}
			pMarginalDistr[i1].MarginalDistrSNP_Y[i2][0] = count;
			// pMarginalDistr[i1].MarginalDistrSNP_Y[i2][0] stores the number in the cell X_i1 = i_2 and Y = 1

			count = 0;
			for (i3 = 0; i3< nlongintctrl; i3++)
			{
				count += bitCount(pgeno[i1*3 + i2].genoctrl[i3]);
			}
			pMarginalDistr[i1].MarginalDistrSNP_Y[i2][1] = count;
			// pMarginalDistr[i1].MarginalDistrSNP_Y[i2][1] stores the number in the cell X_i1 = i_2 and Y = 0

			pMarginalDistr[i1].MarginalDistrSNP[i2] = pMarginalDistr[i1].MarginalDistrSNP_Y[i2][0] + pMarginalDistr[i1].MarginalDistrSNP_Y[i2][1];
        }
    }
}

// GenoDistr[i1*3 + i2] stores the number of cases in the cell X_j1 = i1 & X_j2 = i2
// GenoDistr[9 + i1*3 + i2] stores the number of corresponding controls
void CalculateGenoJointDistr(struct genotype *pgeno, int nsnp, int nLongIntcase, int nLongIntctrl, int *GenoDistr,
	int j1, int j2, struct MarginalDistr *pMarginalDistr)
{
	int i1, i2, i3;

	register int count;
//	uint64 tmp;

	for (i1 = 0; i1<2 ; i1++)
	{
		for (i2 = 0; i2 <2; i2++)
		{
		    count = 0;
			for (i3 = 0; i3< nLongIntcase; i3++)
			{
			    //count += bitCount(pgeno[j1*3 + i1].genocase[i3] & pgeno[j2*3 + i2].genocase[i3]);
				count += popcount(pgeno[j1*3 + i1].genocase[i3] & pgeno[j2*3 + i2].genocase[i3]);
				//GenoDistr[i1*3 + i2] += bitCount(pgeno[j1*3 + i1].genocase[i3] & pgeno[j2*3 + i2].genocase[i3]);
			}

            GenoDistr[i1*3 + i2] = count;
            count = 0;

			for (i3 = 0; i3< nLongIntctrl; i3++)
			{
			    //count += bitCount(pgeno[j1*3 + i1].genoctrl[i3] & pgeno[j2*3 + i2].genoctrl[i3]);
				count += popcount(pgeno[j1*3 + i1].genoctrl[i3] & pgeno[j2*3 + i2].genoctrl[i3]);
				//GenoDistr[9 + i1*3 + i2] += bitCount(pgeno[j1*3 + i1].genoctrl[i3] & pgeno[j2*3 + i2].genoctrl[i3]);
			}
			GenoDistr[9 + i1*3 + i2] = count;
		}

	}
//for case
	GenoDistr[2] = pMarginalDistr[j1].MarginalDistrSNP_Y[0][0] - GenoDistr[0] - GenoDistr[1];
	GenoDistr[5] = pMarginalDistr[j1].MarginalDistrSNP_Y[1][0] - GenoDistr[3] - GenoDistr[4];

	GenoDistr[6] = pMarginalDistr[j2].MarginalDistrSNP_Y[0][0] - GenoDistr[0] - GenoDistr[3];
	GenoDistr[7] = pMarginalDistr[j2].MarginalDistrSNP_Y[1][0] - GenoDistr[1] - GenoDistr[4];

	GenoDistr[8] = pMarginalDistr[j2].MarginalDistrSNP_Y[2][0] - GenoDistr[2] - GenoDistr[5];

//for ctrl
	GenoDistr[11] = pMarginalDistr[j1].MarginalDistrSNP_Y[0][1] - GenoDistr[9] - GenoDistr[10];
	GenoDistr[14] = pMarginalDistr[j1].MarginalDistrSNP_Y[1][1] - GenoDistr[12] - GenoDistr[13];

	GenoDistr[15] = pMarginalDistr[j2].MarginalDistrSNP_Y[0][1] - GenoDistr[9] - GenoDistr[12];
	GenoDistr[16] = pMarginalDistr[j2].MarginalDistrSNP_Y[1][1] - GenoDistr[10] - GenoDistr[13];

	GenoDistr[17] = pMarginalDistr[j2].MarginalDistrSNP_Y[2][1] - GenoDistr[11] - GenoDistr[14];

}

// Calculate the logliner statistics using the joint distribution.
double CalcLLStaWithJointDistr(int *GenoJointDistr, int ncase, int nctrl)
{
	int n_i_j_dot;
	double ls_lb = 0;
	int i, j;
	int n = ncase + nctrl;
	for (i = 0; i<3; i++)
		{
			for(j = 0; j<3; j++)
			{
			    n_i_j_dot = GenoJointDistr[i + 3*j] + GenoJointDistr[i + 3*j + 9];
			    ls_lb = ls_lb + GenoJointDistr[i + 3*j] * log((double)(GenoJointDistr[i + 3*j]) * n / n_i_j_dot / ncase)
					+ GenoJointDistr[i + 3*j + 9] * log((double)(GenoJointDistr[i + 3*j + 9]) * n / n_i_j_dot / nctrl);

			}

		}
		//ls_lb = n * ls_lb;      // It's the loglinear statistic value
		return ls_lb;
}

// Calculate the marginal entropy
void CalculateMarginalEntropy(struct genotype *pgeno, int nsnp, int n, int nlongintcase, int nlongintctrl, double *MarginalEntropySNP, double *MarginalEntropySNP_Y)
{
    int i1, i2, i3;
    int count;
	double tmp, ptmp;
    int GenoMarginalDistr[3][2];

    for (i1 = 0; i1< nsnp; i1++)
    {
        for (i2 = 0; i2<3; i2++)
        {
            count = 0;
            for (i3 = 0; i3< nlongintcase; i3++)
			{
			    count += bitCount(pgeno[i1*3 + i2].genocase[i3]);
			}
			GenoMarginalDistr[i2][0] = count;		// The number of cases in the cell X_i1 = i2

			count = 0;
			for (i3 = 0; i3< nlongintctrl; i3++)
			{
			    count += bitCount(pgeno[i1*3 + i2].genoctrl[i3]);
			}
			GenoMarginalDistr[i2][1] = count;		// The number of controls in the cell X_i1 = i2
        }

		for (i2 = 0; i2<3; i2++)
		{
			tmp = (double) GenoMarginalDistr[i2][0] + GenoMarginalDistr[i2][1];		// The number of samples whose X_i1 = i2
			if ( tmp > 0)
			{
				ptmp = tmp/n;		// The proportion of samples whose X_i1 = i2
				MarginalEntropySNP[i1] += -(ptmp)*log(ptmp);
			}

			if (GenoMarginalDistr[i2][0]>0)
			{
				ptmp = (double) GenoMarginalDistr[i2][0]/n;
				MarginalEntropySNP_Y[i1] += -ptmp*log(ptmp);
			}

			if (GenoMarginalDistr[i2][1]>0)
			{
				ptmp = (double) GenoMarginalDistr[i2][1]/n;
				MarginalEntropySNP_Y[i1] += -ptmp*log(ptmp);
			}

		}
    }

}

void marginAssocDetection(struct genotype *pgeno, int p, int ncase, int n, int nlongintcase, int nlongintctrl) {
	double *MarginalEntropySNP, *MarginalEntropySNP_Y, *MarginalAssociation;
	double MarginalEntropyY;
	double ptmp1 = (double) ncase/n;
	int i;
	FILE *fp;
	MarginalEntropyY =  -ptmp1 *log(ptmp1) - (1-ptmp1) *log(1-ptmp1);

	MarginalEntropySNP = (double *)calloc(p, sizeof(double));
	MarginalEntropySNP_Y = (double *)calloc(p, sizeof(double));
	MarginalAssociation = (double *)calloc(p, sizeof(double));

    CalculateMarginalEntropy(pgeno, p, n, nlongintcase, nlongintctrl, MarginalEntropySNP, MarginalEntropySNP_Y);


    for (i = 0; i<p; i++)
    {
        MarginalAssociation[i] = (-MarginalEntropySNP_Y[i] + MarginalEntropySNP[i] + MarginalEntropyY)*n*2;
    }

    fp = fopen("MarginalAssoc.txt","w");
	if(!fp) {
		printf("failed to open output file MarginalAssoc.txt\n");
		exit(1);
	}
    for (i = 0; i < p; i++)
    {
		// Calculate p value
		/*if(MarginalAssociation[i] < 32 && MarginalAssociation[i] >20 )	{
			double p_val = exp(-0.5 * MarginalAssociation[i]);
			fprintf(fp,"%7d %f %le\n", i,MarginalAssociation[i], p_val);
		}*/
		double p_val = exp(-0.5 * MarginalAssociation[i]);
		fprintf(fp,"%7d %f %le\n", i,MarginalAssociation[i], p_val);
    }
    fclose(fp);

	free(MarginalEntropySNP);
	free(MarginalEntropySNP_Y);
	free(MarginalAssociation);
}

// Calculate the loglinear statistics and store them in a list
// pp: the number of the wanted snp pairs
// SNPpairs: stores the indices of the wanted snp pairs, whose size should be pp*2
// statistics: stores the result, whose size should be pp because we have a statistic value for each pair
// GenoJointDistr: stores the joint distribution tables, whose size should be 9*2*pp
// GenoJointDistr[i*pp + j] represents the j-th cell in the table of the i-th snp pair.
void CalcLLStaWithBits(struct genotype *pgeno, int p, int ncase, int nctrl, int nlongintcase, int nlongintctrl)
{
    int k, j1, j2, i, j;
    int n_i_j_dot;
    int tmpp;
	int n = ncase + nctrl;
    double ls_lb = 0, ls_lb_case = 0, ls_lb_ctrl = 0;
	double ls_lp1 = 0, ls_lp_case1 = 0, ls_lp_ctrl1 = 0;
	double ls_lp2 = 0, ls_lp_case2 = 0, ls_lp_ctrl2 = 0;
    double sta;
	int df;
	double bic_f, bic_p1, bic_p2;

    struct MarginalDistr *pMarginalDistr = (struct MarginalDistr *)malloc(sizeof(struct MarginalDistr) * p);
    int *GenoJointDistr = (int *)calloc(18, sizeof(int));

	char *output_filename = "interaction.out";
	FILE *f_output;
	f_output = fopen(output_filename, "w");
	if(f_output == NULL)
	{
		fprintf(stderr,"can't open output file %s\n", output_filename);
		return;
	}

    // calculate marginal distribution
	CalculateMarginalDistr(pgeno, p, n, nlongintcase, nlongintctrl, pMarginalDistr);
	//free(pMarginalDistr);
	for(j1 = 0; j1 < p; j1 += 2)
	{
		j2 = j1 + 1;
		/*tmpp = j1;
		j1 = j2;
		j2 = tmpp;*/
        // calculate joint distribution for SNP j1 and SNP j2
        CalculateGenoJointDistr(pgeno, p, nlongintcase, nlongintctrl, GenoJointDistr, j1, j2, pMarginalDistr);
        for (i = 0; i<3; i++)
        {
            for(j = 0; j<3; j++)
            {
                n_i_j_dot = GenoJointDistr[i + 3*j] + GenoJointDistr[i + 3*j + 9];
     //           if(GenoJointDistr[i + 3*j] < 1 || GenoJointDistr[i + 3*j + 9] < 1)
     //           {
     //               //printf("n_ij == 0!");
     //               //ls_lb = -1;
					////break;
					//continue;
     //           }
				ls_lb_case = 0;
				ls_lb_ctrl = 0;
				ls_lp_case1 = 0;
				ls_lp_ctrl1 = 0;
				ls_lp_case2 = 0;
				ls_lp_ctrl2 = 0;
				if(GenoJointDistr[i + 3*j] > 0) {
					ls_lb_case = GenoJointDistr[i + 3*j] * log((double)(GenoJointDistr[i + 3*j]) * n / n_i_j_dot / ncase);
					ls_lp_case1 = GenoJointDistr[i + 3*j] * log((double)(GenoJointDistr[i + 3*j]) * pMarginalDistr[j1].MarginalDistrSNP[i] / n_i_j_dot / pMarginalDistr[j1].MarginalDistrSNP_Y[i][0]);
					ls_lp_case2 = GenoJointDistr[i + 3*j] * log((double)(GenoJointDistr[i + 3*j]) * pMarginalDistr[j2].MarginalDistrSNP[i] / n_i_j_dot / pMarginalDistr[j2].MarginalDistrSNP_Y[i][0]);
				}
				if(GenoJointDistr[i + 3*j + 9] > 0) {
					ls_lb_ctrl = GenoJointDistr[i + 3*j + 9] * log((double)(GenoJointDistr[i + 3*j + 9]) * n / n_i_j_dot / nctrl);
					ls_lp_ctrl1 = GenoJointDistr[i + 3*j + 9] * log((double)(GenoJointDistr[i + 3*j + 9]) * pMarginalDistr[j1].MarginalDistrSNP[i] / n_i_j_dot / pMarginalDistr[j1].MarginalDistrSNP_Y[i][1]);
					ls_lp_ctrl2 = GenoJointDistr[i + 3*j + 9] * log((double)(GenoJointDistr[i + 3*j + 9]) * pMarginalDistr[j2].MarginalDistrSNP[i] / n_i_j_dot / pMarginalDistr[j2].MarginalDistrSNP_Y[i][1]);
				}
				ls_lb += ls_lb_case + ls_lb_ctrl;
				ls_lp1 += ls_lp_case1 + ls_lp_ctrl1;
				ls_lp2 += ls_lp_case2 + ls_lp_ctrl2;
                /*ls_lb += GenoJointDistr[i + 3*j] * log((double)(GenoJointDistr[i + 3*j]) * n / n_i_j_dot / ncase)
                    + GenoJointDistr[i + 3*j + 9] * log((double)(GenoJointDistr[i + 3*j + 9]) * n / n_i_j_dot / nctrl);*/

            }
			/*if(ls_lb == -1) {
				break;
			}*/

        }
        //ls_lb = n * ls_lb;      // It's the loglinear statistic value

        // Do Chi-square test with df = 8
        // Maybe we donnot need to do the test, just count the rank
		bic_f = ls_lb + 8*log(n);
		bic_p1 = ls_lp1 + 6 * log(n);
		bic_p2 = ls_lp2 + 6 * log(n);
		if(bic_f <= bic_p1 && bic_f <= bic_p2) {
			df = 8;
			sta = ls_lb;
		}
		else if(bic_p1 <= bic_f && bic_p1 <= bic_p2) {
			sta = ls_lp1;
			df = 6;
		}
		else if(bic_p2 <= bic_f && bic_p2 <= bic_p1) {
			df = 6;
			sta = ls_lp2;
		}
		sta = sta * 2;
        ls_lb = 0;
		ls_lp1 = 0;
		ls_lp2 = 0;
			
		// write to output file
		if(sta > 72.43) {
			// p_value = 1.61e-12, sta = 72.43
			// p_value = 10^-6, sta = 42.7
			// p_value = 10^-13, sta = 78.45
			//fprintf(f_output, "%d\t%d\t%f\n", j1, j2, sta);
		}
		fprintf(f_output, "%d\t%d\t%f\t%d\n", j1, j2, sta, df);
		/*tmpp = j1;
		j1 = j2;
		j2 = tmpp;*/
    }
	fclose(f_output);
	free(pMarginalDistr);
	free(GenoJointDistr);
	return;
}

void PermuteOnContingencyTable(int *data_matrix, int *labels, int *perm_labels, int *SNPpairs, int p, int ncase, int nctrl, int pp,
	int *GenoJointDistr, double *sta, int *perm_cnt)
{
	double perm_sta;
    int i, j, k;
    int tmp;
	int n = ncase + nctrl;
	int p1, p2;
	int g1, g2;
	int delta;

	memcpy(perm_labels, labels, n * sizeof(int));

	for(i = 0; i < n; i++)
	{
	    j = random(n - i) + i;
	    if(perm_labels[i] != perm_labels[j])
		{
			// swap the i-th and the j-th label
			tmp = perm_labels[i];
			perm_labels[i] = perm_labels[j];
			perm_labels[j] = tmp;
		}
	}
	// Now we get the permuted labels in perm_labels[]

    for(i = 0; i < n; i++)
    {
		if(perm_labels[i] != labels[i])
		{
			// modify the contingency table
			for(k = 0; k < pp; k++)
			{
				p1 = SNPpairs[2*k];
				p2 = SNPpairs[2*k + 1];
				// X_p1 = data_matrix[i * p + p1] && X_p2 = data_matrix[i * p + p2]
				// <---->exchange with
				// X_p1 = data_matrix[j * p + p1] && X_p2 = data_matrix[j * p + p2]
				g1 = data_matrix[i * p + p1];
				g2 = data_matrix[i * p + p2];
				delta = perm_labels[i] == 1 ? 1 : -1;
				GenoJointDistr[18 * k + g1 * 3 + g2] += delta;
				GenoJointDistr[18 * k + 9 + g1 * 3 + g2] -= delta;
			}
		}
    }
	for(k = 0; k < pp; k++)
	{
		// The permutation finished. Start to calculate the new statistics.
		perm_sta = CalcLLStaWithJointDistr(GenoJointDistr + 18 * k, ncase, nctrl);
		if(perm_sta > sta[k])
		{
			perm_cnt[k]++;
		}
	}

}

int main()
{
	time_t st, ed;

	struct genotype *pgeno;
	int *SNPpairs; // the indices of the wanted snp pairs.
	int *labels, *perm_labels, *tp;

	/* Declare variable */
	int *DataSize;
	int ndataset;

    int n, p;  //n--number of samples; p number of varibles
	int ncase, nctrl, nlongintcase, nlongintctrl;
	int icase, ictrl;

	int flag, i, ii, j, k;
	int tmp;
	int cnt = 0;

	int *GenoJointDistr;
	double *original_sta;
	//double *perm_sta;
	int *perm_cnt;

	// used for association detection allowing for interaction

	int LengthLongType=64;
	uint64 mask1 = 0x0000000000000001;

	// load data
	FILE *fp_pairs;
	FILE *fp, *fp_i;
	char filename_i[100];
	char filename[100] = "filenamelist.txt";
	char filename_pairs[100] = "pairs.txt";
	fp = fopen(filename,"r");
	if(fp == NULL)
	{
		fprintf(stderr,"can't open input file %s\n",filename);
		return 1;
	}

	// precompute the wordbits (a global variable)
    for (i = 0; i<65536; i++)
    {
        wordbits[i] = bitCount(i);
        //printf("%d\n",wordbits[i]);
    }
	printf("start loading ...\n");

	ndataset = GetDataSize(filename, &DataSize);

	n = DataSize[0];
	p = 0;
	printf("n = %d\n", n);
	for (i = 0; i<ndataset ; i++)
	{
		p += DataSize[ndataset*1 + i];
		printf("DataSize %d-th file: p[%d] = %d \n", i+1, i+1, DataSize[ndataset*1 + i]);
	}
	printf("p = %d\n",p);
	labels = (int *)calloc(n, sizeof(int));
	perm_labels = (int *)calloc(n, sizeof(int));

	// get ncase and nctrl
	i = 0;
	j = 0;

	ncase = 0;
	nctrl = 0;

	rewind(fp);

	// only use the first file to get ncase and nctrl
	fscanf(fp, "%s\n", &filename_i);
	printf("%s\n", filename_i);
	fp_i = fopen(filename_i, "r");

	while(!feof(fp_i)) {

		/* loop through and store the numbers into the array */
		if(j==0)
		{
			//j = 0 means read ind class label y
			fscanf(fp_i, "%d", &tmp);

			if (tmp)
			{
				// tmp=1 means case
				ncase++;

			}
			else
			{
				nctrl ++;

			}
			j++;
		}
		else
		{
			fscanf(fp_i, "%d", &tmp);
			j++; //column index
			if (j==(DataSize[ndataset]+1)) // DataSize[ndataset] is the nsnp in the first dataset
			{
				j=0;
				i++; // row index
			}

		}

		if (i>=n)
		{
			break;
		}
	}


	printf("total sample: %d (ncase = %d; nctrl = %d).\n", n, (int)ncase, (int)nctrl);

	nlongintcase = ceil( ((double) ncase)/LengthLongType);
	nlongintctrl = ceil( ((double) nctrl)/LengthLongType);
	printf("nLongIntcase = %d; nLongIntctrl = %d.\n", nlongintcase, nlongintctrl);

	//calloc memory for bit representation
	pgeno = (struct genotype *)malloc(sizeof(struct genotype) * p * 3);// p SNPs, each contains 3 genotypes
	for (j = 0; j< 3*p ; j++)
	{
		(pgeno+j)->genocase = (uint64 *)calloc( nlongintcase, sizeof(uint64));
		(pgeno+j)->genoctrl = (uint64 *)calloc( nlongintctrl, sizeof(uint64));
	}

	//load data to bit representation


	rewind(fp);

	time(&st);
	j = 0; // column index
	ii = 0; // file index
	k = 0;
	while(!feof(fp)) {
		ii ++;
		fscanf(fp, "%s\n", &filename_i);

		fp_i = fopen(filename_i, "r");
		if(fp_i == NULL)
		{
			fprintf(stderr,"can't open input file %s\n", filename_i);
			exit(1);
		}

		i = 0; //row index
		icase = -1;
		ictrl = -1;

		printf("Loading data in file %d: %s\n", ii, filename_i);
		while(!feof(fp_i)) {
			/* loop through and store the numbers into the array */

			if(j==0)
			{
				//j = 0 means read class label y
				fscanf(fp_i, "%d", &tmp);
				// stores the labels for permutation
				labels[i] = tmp;
				if (tmp)
				{
					// tmp=1 means case
					icase ++;
					flag = 1;
				}
				else
				{
					ictrl ++;
					flag = 0;
				}
				j++;
			}
			else
			{
				fscanf(fp_i, "%d", &tmp);
				// stores the snp data for permutation
				if (flag)
				{
					pgeno[(j+k-1)*3 + tmp].genocase[icase/LengthLongType] |= (mask1 << (icase%LengthLongType));
				}
				else
				{

					pgeno[(j+k-1)*3 + tmp].genoctrl[ictrl/LengthLongType] |= (mask1 << (ictrl%LengthLongType));
				}

				j++; //column index
				if (j==(DataSize[ndataset + ii-1]+1))
				{
					j=0;
					i++; // row index
				}

			}

			if (i>=n)
			{
				break;
			}
		}

		fclose(fp_i);
		k += DataSize[ndataset + ii-1];
	}

	fclose(fp);

	//The old solution is to use a pair index file to specify the pairs. 
	//The new design is that in the genotype data file, the SNPs are grouped by pairs. 
	//i.e. SNP0 and SNP1 is an interacted pair. So the total number of SNPs should be even. 

	// Calculate the log linear statistics for the original data.
	CalcLLStaWithBits(pgeno, p, ncase, nctrl, nlongintcase, nlongintctrl);
	marginAssocDetection(pgeno, p, ncase, n, nlongintcase, nlongintctrl);

	free(DataSize);
	for (j = 0; j< 3*p ; j++)
	{
		free(pgeno[j].genocase);
		free(pgeno[j].genoctrl);
	}

	free(pgeno);
	free(labels);
	free(perm_labels);
	//free(GenoJointDistr);
	//free(original_sta);

	return 1;
}
